      subroutine initxsec

      use param
      use allocs
      
#include <global.h>
	include 'ntbytes.h'
      include 'geom.h'
      include 'geomf.h'
      include 'xsec.h'
      include 'cmfd.h'
      include 'cmfdp1.h'

      real(NBF) coef

      coef=3./7

      allocate (xssm(tnode,ng))
      allocate (xsmf(fnode,ng))

      do im=1,ng
        do i=1,tnode
          allocate(xssm(i,im)%from(ng))
        enddo
        do j=1,fnode
          allocate(xsmf(j,im)%from(ng))
        enddo
      enddo

!     transfer group/composition xsec to region/node xsec
      do im=1,ng
        do j=1,tnode
          icomp=comp(j)
          xsd(1,j,im)=diff(icomp,im)  ! D0
          xsd(2,j,im)=coef*rsigt(icomp,im)    ! D2 
          xst(j,im)=sigt(icomp,im)
          xsr(j,im)=sigr(icomp,im)
          xsnf(j,im)=signf(icomp,im)
          xchi(j,im)=chi(icomp,im)
          xssm(j,im)%from(1:ng)=sigsm(icomp,im)%from(1:ng)
          rxst(j,im)=1/xst(j,im)
          rxsr(j,im)=1/xsr(j,im)
        enddo
!
        do j=1,fnode
          icf=compf(j)
          xsdf(1,j,im)=diff(icf,im)  ! D0
          xsdf(2,j,im)=coef*rsigt(icf,im)    ! D2 
          xstf(j,im)=sigt(icf,im)
          xsrf(j,im)=sigr(icf,im)
          xsff(j,im)=signf(icf,im)
          xchif(j,im)=chi(icf,im)
          xsmf(j,im)%from(1:ng)=sigsm(icf,im)%from(1:ng)
          rxstf(j,im)=1/xstf(j,im)
          rxsrf(j,im)=1/xsrf(j,im)
        enddo
      enddo

!     define beta
!     left(bottom) boundary
      if (bcb.eq.0) then
        bcm_p1(1)=0
        bcm_p3(1,:)=bcref
      elseif (bcb.eq.1) then
        bcm_p1(1)=bcin0_p1(1)
        bcm_p3(1,:)=bcin0_p3
      endif
!     right(upper) boundary
      if (bcu.eq.0) then
        bcm_p1(2)=0
        bcm_p3(2,:)=bcref
      elseif (bcu.eq.1) then
        bcm_p1(2)=bcin0_p1(2)
        bcm_p3(2,:)=bcin0_p3
      endif

      do im=1,ng
        betax(:,0,im)=2*beta(:,0,im)
        do j=1,tnode
          beta(:,j,im)=xsd(:,j,im)*rhz(j)
          betax(:,j,im)=2*beta(:,j,im)
        enddo
        betax(:,j,im)=2*beta(:,j,im)
      enddo

!     finemesh beta
      do im=1,ng
        do j=1,fnode
          betaf(:,j,im)=xsdf(:,j,im)*rhzf(j)
        enddo
      enddo

      end subroutine
